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The ground state energies of infinite half-filled Hubbard- 
Peierls chains are investigated combining incremental expan- 
sion with exact diagonalization of finite chain segments. The 
ground state energy of equidistant infinite Hubbard ( Heisen- 
berg) chains is calculated with a relative error of less than 
3 • 10"^ for all values of U using diagonalizations of 12-site 
(20-site) chain segments. For dimerized chains the dimer- 
ization order parameter d as a function of the onsite repul- 
sion interaction U has a maximum at nonzero values of [7, if 
the electron-phonon coupling g is lower than a critical value 
Qc- The critical value gc is found with high accuracy to be 
gc = 0.69. For smaller values of g the position of the maxi- 
mum of d{U) is approximately 3t, and rapidly tends to zero 
as g approaches g^ from below. We show how our method can 
be applied to calculate breathers for the problem of phonon 
dynamics in Hubbard-Peierls systems. 



The effect of correlations on the Peierls transition has 
been one of challenging problems in the theory of quasi- 
one-dimensional compounds. One of the most impor- 
tant theoretical treatments of the Peierls transition goes 
back to the solution of the exactly solvable model of 
noninteracting fermions proposed by Su, Schrieffer and 
Heeger (SSH) ^ Although being successful in explain- 
ing a number of properties of real quasi-one-dimensional 
systems, the SSH model is in a clear disagreement with 
such experimental results as the emergence of negative 
spin magnetization densities for neutral solitons . One 
is faced with a necessity to treat the Coulomb interac- 
tion in the electron subsystem. This interaction should 
be accounted for by including a positive Hubbard onsite 
interaction term in the SSH model. We refer to this ex- 
tended model as to Peierls-Hubbard model (PHM). Due 
to strong one-dimensional quantum fluctuations a mean 
field theory calculation of PHM gives qualitatively wrong 
results, predicting a constant dimerization for small U 
and the abrupt disappearance of a bond order wave state 
upon increasing U above a certain threshold at half-filling 
1^. Including many-body effects it has been shown by 
many authors ( Q and citations therein) that the dimer- 
ization d first increases up to a maximum and then de- 
creases with further increase of U. 

It is very difficult to perform an accurate exact diag- 
onalization investigation of the Peierls transition in the 
correlated regime . In the framework of a standard exact 
diagonalization approach the required cluster sizes are 
found to be far outreaching the capabilities of modern 



computer systems. Calculations which were performed 
using available cluster sizes are drastically depending on 
the boundary conditions (see for instance ||^ ) and the 
final conclusions had to be based on extrapolations. The 
basic questions about the value of the behaviour of 
the system near the critical point, the position and the 
value of the dimerization maximum Umax as a function 
of the electron phonon coupling remained unanswered. 
The lack of accurate numerical results was making it hard 
to identify the values of model parameters for real sys- 
tems. We show that by combining an incremental expan- 
sion technique (lET) with numerically exact diagonaliza- 
tions, one can overcome the abovementioned difficulties 
and perform a reliable numerical calculation of correlated 
one-dimensional Peierls systems in both strong and weak 
correlation regimes. 

The quantum chemical method of increments found 
recently a wide region of application in condensed mat- 
ter (see Q and references therein). The lET starts 
with splitting a given Hamiltonian operator H into an 
unperturbed part Hq and a number of perturbations 
Hi+ H2 + H = HQ + J2Hi. The hierarchy of in- 
crements is defined in the following way. The first order 
increment ij.^^ is given by taking the ground state energy 
£k of the Hamiltonian Hq + Hk and subtracting from it 
the ground state energy £0 of Hq, /^^^ = £k — So- Phe- 
nomenologically this increment represents the action of 
the perturbation Hk separately from all other perturba- 
tions Hi. The total change of the unperturbed ground 
state energy in first incremental order is then given by 
The second order increment I^^ is defined by 
taking the ground state energy £k,i of the Hamiltonian 
Ho+Hk+Hi and subtracting from it Sq o,f^d the first order 



increments I^'^ and ^^^"^ 



^k,l 



r(i)_f^(i). 



-fn. The 



increments / represent the difference between the com- 
bined action of a pair Hi, Hk and the sum of the uncorre- 
lated actions of both perturbations. In a similar manner 
higher order increments are defined. The change of the 
unperturbed ground state energy of the full system is ex- 
actly given by the sum over all increments (all orders!). 
Since the increments are usually calculated numerically, 
the incremental expansion can be performed up to some 
given order. This expansion is nonperturbative since the 
increments are not related to some small parameter of 
a perturbation theory. The idea of the incremental ex- 
pansion is similar to Faddeev's treatment of the 3-body 
problem where the unknown 3-body scattering matrix is 
expressed through the exactly known 2-body scattering 
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matrices. The discussion of the interrelation of incre- 
ments and Faddeev equations and also the derivation of 
the incremental expansion by a resummation of the per- 
turbation theory is given in ( ||^). 

Now we apply the outlined ideas to the PHM Hamil- 
tonian. This Hamiltonian is given by a sum of electronic 
and lattice parts H = H^i -\- Hiat- The electronic part in 
fcrmionic second quantization form is given by 

-ffe; = ^tj(4o.Cj+i,o- -t- /i.e.) + [/^(n^jn^^x) . (1) 

Here U is an onsite Hubbard repulsion matrix element 
and ti is the hopping matrix element between the i-ih 
and i + 1-th sites. We consider the case of one electron 
per site (half filling). In the harmonic approximation 
the lattice part reads as Hiat = 1/2-?^ X^i ■ Here Vi is a 
bond-length change (see e.g. Q) and K is the spring con- 
stant. The electron-lattice interaction is assumed to be 
of the form ti = —{t — ^Vi). The strength of the electron- 
phonon interaction is measured by the dimensionless cou- 
pling g = 7/%/ Kt. Solving the PHM amounts to finding 
a minimum of the full energy of the system considered 
as a functional of the bond length changes Vi. A remark- 
able proof of Lieb and Nachtergaele Q tells us that the 
minimum configuration has to be a dimerized state with 
alternating bond lengths Vi — (— l)*wo. In the following 
we will use the dimensionless dimerization d = voy^K/t 
(see ref. 1,1). 

Let us now formulate the incremental expansion of the 
PHM. A dimerised state represents a sequence of alter- 
nating weak and strong bonds formed by a modulation 
of the transfer integral ti. It is natural to cut all the 
weak bonds and to consider the remaining set of nonin- 
teracting 2-site dimers as an unperturbed Hamiltonian 
Hq. The weak bonds are considered as a perturbation. 
The unperturbed Hamiltonian Hq is written as 

00 00 

k— — 00 i — — OQ,(T 

(2) 

The ground state of the Hamiltonian operator Hq is 
known exactly and is a nondegenerate spin-singlet state 
S — Q formed by a set of noninteracting dimers hav- 
ing two electrons per dimer. The PHM hamiltonian H 
is a sum of Hq and a number of perturbations, formed 
by the (weak) bonds linking neighbouring dimers H = 

Ho + J2k^k, Vk = J2a'^2k-l{c2k-l,a^2k.cr + h.c). The 

incremental expansion is generated in the following way. 
The first order increment corresponds to a bond inserted 
between two neighbouring dimers. Due to general prin- 
ciples outlined above it reads as: /^^^ = E^ — 2E2 ; where 
E2n, n = 1,2, ... denotes the ground state energy of a 2n 
sites segment cut out of the infinite chain. The second 
order increment is defined for a triple of neighbouring 



dimers and follows from inserting two bonds into Hq. To 
find it one needs to subtract from the energy of three 
linked dimers the increments corresponding to 2 pairs of 
dimers and 3 single dimers in it, hence the expression 
reads as: J^^) ^ Ee- 2/(i) - 3E2 =^ Ee~ 2Ei + E2. Note 
that in the incremental expansion only connected clus- 
ters of dimers yield nonzero increments, since the energy 
of a disconnected cluster is just the sum of the energies 
of its parts. The expressions for higher order increments 
are found in similar fashion. Due to the choice of Hq 
(see (|^) ) the increments do not depend on site indices. 
One proves by induction that the expression for the n-th 
order increment, n > 2 reads as 

/(") = E2n+2 — 2i?2n + E2n-2 ■ (3) 

In order to find the ground state energy of the infinite 
system one needs to count the order of increments of 
each order per dimer. In the infinite 1-d lattice which 
we consider, there exists exactly one increment of each 
order per dimer (one makes a one-to-one correspondence 
between dimers and increments of each order, assigning 
each increment to its left most dimer). Therefore, taking 
into account that there are two sites per each dimer, the 
value of the ground state energy of the infinite lattice per 
site is written as f = l/2(£;2 + E^=i-^^"-')- Here N is 
the number of increments taken into account. 

The feature that the number of increments of any given 
order per site is constant is an exclusive property of one- 
dimensional lattices. In higher dimensions the per site 
number of increments of a given order grows rapidly with 
the order of the increment. This special property of one 
dimension leads to the following result: 

£{N) = ^{E2N+2-E2n) ■ (4) 

Formula (^) is quite remarkable, since the calculation of 
ground state properties amounts to the exact diagonal- 
ization of two open chains whose length differs by two. 
Note that expressions of the type (^) were intuitively used 
before in quantum chemical calculations (see for instance 
§)■ 

To check our method we first performed a calculation 
of the ground state energy of an equidistant Hubbard in- 
finite chain at half filling, where the solution is known 
exactly [|l^. The equidistant case is the worst case for 
the method described above, since all the bonds have 
the same strength. The per site value of the ground 
state energy £ was calculated using the formula with 
the incremental order N — 1,2,3,4,5. The calculation 
was performed using a Lanczos algorithm. The results 
for iV = 5 are shown in the Fig. 1, where the exact 
£{U) dependence for t = 1 (solid line) is plotted against 
the results of (jj) (open circles). The relative error Re 
decays algebraically with growing order of increments 
Re = A{U)[2N + 2]-"^^'). The exponent viU) and the 
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prefactor A{U) are plotted in the inset of Fig.l. We find 
i^{U) > 2 for all values of U. Note that A{0)/A{oo) w 2 
which implies that our results converge faster for large 
U. Note also that the errors are very small - for N = 6 
typically below 0.1%. This has to be compared with a 
recent Density Matrix Renormalization Group (DMRG) 
calculation of the same system [Q, where system sizes 
up to 122 and extrapolations had to used to achieve com- 
parable precision. 

Next we show the results of calculations of the dimen- 
sionless dimerization in the PHM as a function of U and 
g (Fig. 2). For U — the value of d{g) is known exactly 
(see filled symbols at ?7 = in Fig. 2). An analysis 
of the relative error Rd of determining d with the help 
of (Q) yields exponential convergence Rd ~ 6"'^'^^)'^^+^^ 
The dependence of X{g) is shown in the inset of Fig. 2. A 
crossover is detected around g — 0.4 with A being sup- 
pressed to rather small values for g < 0.4. That implies 
that for small values of U the lET method using exact 
diagonalizations is confined to values of the coupling con- 
stant g > 0.4 if sufficient precision is requested. In Fig. 2 
we present the dependence of d{U) for g = 0.5, 0.6 and 
0.7 (open symbols). For g = 0.5, 0.6 the dimerization d 
first increases with U, and then decreases after reaching 
a maximum. For g — 0.7 is a monotonously decreasing 
function of U. Therefore the system has a qualitatively 
different behaviour for weak and strong couplings g, as 
predicted by the GA theory. On the other hand, our 
results well agree with the extrapolated values of d ob- 
tained within the solitonic approach |Q. 

We performed more calculations of d{U) to obtain the 
dependence of the position of the maximum Umax on the 
coupfing g (see Fig. 3). Especially we find Umax{g) to be 
a monotonously decreasing function with Umax = at a 
critical coupling gc- Since our method yields very small 
errors for g > 0.6 we can estimate the critical coupling 
g where Umax = with high accuracy Rd < 10"'^. It is 
found gc = 0.69. The GA prediction gc > 0.74 overesti- 
mates this result slightly. The GA gives the position of 
the dimerization maximum as Umax = 4i for g far be- 
low gc- Our numerical calculation gives Umax — St. The 
small U behaviour of d is found to be d ^ U^, which 
is consistent with the GA approach. Furthermore the 
GA approach predicts that d is an analytic function of 
U^. Then it follows that Umax{g) close to gc varies as 
K{gc — gY^"^, which is what we find in Fig. 3 (k ~ 8.25 

For large values of U the system is equivalent to the 
Heisenberg spin-exchange model, with Ji given by Ji = 
Atf/U . We have calculated the spin-Peierls transition in 
this system separately, using the formula (^ and N — 
The results for the dimerization are plotted in Fig. 2 
(filled symbols). Note that the dimerization d{U) for 
the Hubbard and Heisenberg chains converge for large U 
which supports the correctness of our calculations. For 
the spin Peierls transition Inagaki and Fukuyama p^ ] 
found an asymptotic formula 
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(5) 



where D is a constant which was assumed to be of the 
order of 1/2. Our results confirm this choice. 

The above results show that the lET can be a key 
method for numerical studying the static properties of 
one-dimensional Peierls systems. Our recent calculations 
show that it can be equally good applied to spin-Peierls 
systems with frustration We believe also that this 
method could be applied to higher dimensional systems, 
namely to the Peierls transition in the two-dimensional 
Hubbard model. The cancellation of the lower order en- 
ergies does not take place in the two-dimensional Hub- 
bard model so the analogue of the formula (^ contains 
energies of the clusters of all sizes. 

To further underline the applicability of the lET, we 
consider dynamical properties of finite Hubbard-Peierls 
systems. The dynamics of classical degrees of freedom 
(phonons) Qi interacting on a lattice generically allows 
for time-periodic and spatially localized solutions namely 
discrete breathers (DB) if the equations of motion incor- 
porate nonlinear terms (see and citations therein). 
These DB solutions can be localized on as few as three 
neighbouring sites. If the electron-phonon coupling is 
taken into account, and the Born-Oppenheimer approx- 
imation is used, the electronic subsystem generates an 
additional potential for the classical phonon degrees of 
freedom. To numerically find again DB solutions, one 
needs the electronic energy £{{Qi}) as a function of the 
phonon degrees of freedom. For a lattice with L sites 
this amounts to calculating the ground state energy of 
the electronic system L times on each time step in order 
to find the gradient of Hd ■ Precalculating the function 
Hci on a grid is also impossible since it is a function of 
prohibitively many variables. In the static dimerization 
case, where it is known that the target state is a bond 
conjugate state this problem is avoided since there is only 
one variable d. 

Again the lET helps to overcome this problem. Con- 
sider a finite chain with periodic boundary conditions. 
The first order increment I^^^{x), which does not depend 
on the site index, is obtained by fixing all the Qi ~ Q ex- 
cept one with Qi = x, and calculating the change of the 
electronic energy as a function of x, I^^^x) = £{x)-£{0). 
The second order increment is obtained by fixing all 
phonon variables Qi = except Qi = x and Qm — y- 
Then the energy of the electronic system will depend on 
k — I — m , x and y. The second order increment reads 
8. as /f ^ {x, y) = £{x, y) - I^^'> (x) - /(^^ (y) - £{0). In the 
same manner higher order increments are found. Our 
calculations (see |l^] for a detailed discussion) show that 
taking into account increments of first and second order 
is enough to calculate the ground state energy of a 14- 
site Hubbard chain with periodic boundary conditions for 
an arbitrary configuration of {Qi}. The relative error is 
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less than 10~^. The increments are calculated on a two- 
dimensional grid to generate smooth functions. With the 
help of these functions the lattice dynamics can be cal- 
culated using ordinary molecular dynamics techniques. 

In this paper we combined the lET with an exact di- 
agonalization method. It is known that the DMRG is 
especially accurate when applied to large but finite open 
chains, where one can achieve higher and higher accu- 
racy by iteratively repeating the DMRG procedure . 
Taking this into account we think that the combination 
of the lET with the DMRG technique can significantly 
improve calculations. 

We thank Prof. Peter Fulde for fruitful discussions and 
continuous support. 
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Fig.l 

Ground state energy of the equidistant Hubbard model 
at half filling. Solid line - exact result , open circles - 
lET resuh for iV = 5 (see text). 

Inset: J7-dependence of prefactor A and exponent v of 
the found functional dependence of the relative error on 
(27V + 2) (see text). 

Fig.2 

Dimerization versus U for different values of g. Open 
symbols - lET results for iV = 5, filled symbols for U ~ 
- exact results, filled symbols for [/ > 5 - results for 
Heisenberg chains with lET and N = 8. Solid lines are 
guides to the eye. 

Inset: (^-dependence of the exponent A of the found func- 
tional dependence of the relative error on {2N + 2) for 
U = (see text). 



Fig.3 

Dependence of Umax on the coupling g for the Hubbard- 
Peierls chains with lET and N = 5. Circles - numerical 
results; fine - fit (see text). 
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